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We study a problem of non-adiabatic superfluid dynamics of spin-orbit coupled neutral fermions in two 
spatial dimensions. We focus on the two cases when the out-of-equilibrium conditions are initiated either by a 
sudden change of the pairing strength or the population imbalance. For the case of zero population imbalance 
and within the mean-field approximation, the non-adiabatic evolution of the pairing amplitude in a collisionless 
regime can be found exactly by employing the method of Lax vector construction. Our main finding is that the 
presence of the spin-orbit coupling significantly reduces the region in the parameter space where a steady state 
with periodically oscillating pairing amplitude is realized. For the collisionless dynamics initiated by a sudden 
disappearance of the population imbalance we obtain an exact expression for the steady state pairing amplitude. 
In the general case of quenches to a state with finite population imbalance we show that there is a region in 
the steady state phase diagram where at long times the pairing amplitude dynamics is governed by the reduced 
number of the equations of motion in full analogy with exactly integrable case. 

PACS numbers: 05.30.Fk, 32.80.-t, 74.25.Gz 


I. INTRODUCTION 

Starting with the seminal paper by Gor’kov and RashbaP 
there has been a remarkable resurgence of interest in the phys¬ 
ical properties of the spin-orbit coupled superfluids and super¬ 
conductors in the past decade. 2 10 This interest is largerly mo¬ 
tivated by theoretical discovery of topological insulators and 
topological superconductors in which spin-orbit coupling of¬ 
ten plays a crucial role by giving rise to the existence of robust 
conducting states at a system’s boundaries on a background of 
a gapped single particle spectrum in a bulk. 1 3GH In addition, 
the recent discovery of the superconductivity at the interface 
in the oxide-based heterostructures 19 21 where the inversion 
symmetry is naturally broken served as an additional motiva¬ 
tion for studying both conventional and unconventional super¬ 
conductivity in spin-orbit coupled systems. 22 

Of special interest are the physical properties of topologi¬ 
cal insulators and superconductors under external influences 
which drive these systems far-from-equilibrium. In particular, 
the concept of the Floquet topological insulators have been 
recently developed in the context of various systems exter¬ 
nal periodic driving, which leads to an inversion of the bands 
with different parity giving rise to metallic edge states. 23 ® 
Furthermore, several groups have generalized the idea of 
Floquet topological insulators to Floquet topological s-wave 
superconductors ®K] Most recently, it has been shown that 
topological Floquet superfluidity can be realized in systems 
where the periodic driving is self-generated in the process of 
the collisionless dynamics. 32 ® 

However, certain aspects of the pairing dynamics in the col¬ 
lisionless regime for the spin-orbit coupled systems have not 
been addressed yet. The aim of this paper is to close the 
remaining gaps in the studies of this problem. Specifically, 
using both exact integrability and numerical analysis we in¬ 
vestigate how the presence of the spin-orbit coupling affects 
the behavior of the pairing amplitude at long times. We con¬ 
sider the standard protocol of inducing far from equilibrium 
coherent dynamics in fermionic condensates by fast switch 


of one of the system’s parameters. In our model we allow 
for non-zero out-of-plane Zeeman field hz which gives rise 
to the population imbalance between the fermionic atoms in 
two hyperfine states. Here we discuss two cases: changes in 
the detuning frequency of the Feshbach resonance and in the 
population imbalance. 

There are three relevant time scales in the problem: the first 
time is the perturbation time scale T quen ch which we take to 
be instantaneous; the second time scale is governed by the 
dynamics of the Cooper pairs, ta while the third time scale, 
r e , accounts for the relaxation due to two-particle collisions. 
In what follows, we consider the limit r e —>■ oo and analyze 
the dynamics of the pairing amplitude at long times t > r A . 
Importantly, we will also neglect the possibility for the pair¬ 
ing amplitude to become spatially inhomogeneous, which is 
equivalent to an assumption of having a system with a size 
much smaller than the superfluid coherence length. 

Within the mean-field theory for the reduced BCS model in 
the weak coupling limit, three types of steady states have been 
found for the quenches of the pairing strength and provided 
the system is initially in its ground state®® (Regime I) gap¬ 
less steady state with zero pairing amplitude A(t —> oo) = 0; 
(Regime II) steady state with the constant pairing amplitude, 
A(i -A oo) = Aoo; and (Regime III) steady state described by 
the undamped periodic oscillations of the pairing amplitude. 
Interestingly, there are no qualitative changes in the steady 
state phase diagram for the quenches across the s-wave Fes- 
chbach resonanc^l as we fl as f or the two-dimensional chiral 
superfluids. 32 33 In principle, other steady states, such as the 
one in which pairing amplitude is a multiperiod function of 
time, can also be realized. 37 ® However, realization of these 
states requires that the system is initially in an excited state. 

Perhaps the most surprising result from the earlier studies of 
the non-adiabatic pairing problem is the discovery of a steady 
state with the periodically oscillating amplitude whose analyt¬ 
ical expression is given by the Jacobi elliptic function. 35 37 41 42 
Thus, the main thrust of the present work is on one hand to in¬ 
vestigate the fate of that steady state for a condensate with 
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equal populations and non-zero spin-orbit coupling. On the 
other hand, we will also investigate whether in the model with 
the population imbalance a system allows the realization of 
that steady state, i.e. the pairing amplitude is still expressed 
in terms of the Jacobi elliptic function, even though non-zero 
population imbalance precludes the full analytical description. 

Let us briefly summarize our results. In the first part of the 
paper we analyze the effect of the spin-orbit coupling on the 
steady state phase diagram. We find that the steady state III 
is realized in much narrower region of the phase diagram. In 
particular, we find that the size of the Region III is inverse 
proportional to the strength of the spin-orbit coupling. Quali¬ 
tatively, this effect is due to the lifting of the Kramers degen¬ 
eracy by the spin-orbit coupling. Since the total pairing am¬ 
plitude is determined by the pairing in two chiral bands and 
the collective collisionless dynamics is reduced to a motion of 
two effective variables, large spin-orbit coupling effectively 
hinders the appearance of the steady state with periodically 
oscillating amplitude. 

The remaining part of our discussion concerns the nature of 
the steady state for the quenches in the population imbalance. 
This problem has been recently studied by Y. Dong et al. Il34l 
by solving the Bogoliubov-de Gennes equations numerically. 
Here we show that for the quenches to the state with equal 
atomic populations, the superfluid dynamics for the pairing 
amplitude can, in fact, be found exactly. Specifically, we ob¬ 
tain an exact expression for the steady state pairing amplitude 
and analyze the steady state phase diagram as a function of the 
population imbalance in the initial state. Our results for this 
part are generally in agreement with those reported in Ref. 
l34l . Then, we continue with the discussion for the quenches 
to a state with finite population imbalance. For this part we 
had to resort to the numerical analysis of the equations of mo¬ 
tion. Our main finding is that when the finite value of the pop¬ 
ulation imbalance exceeds some critical value, we observe the 
dynamical reduction in the number of quantities describing 
the system’s dynamics. In other words, the order parameter 
dynamics is described by the same equations of motion as in 
integrable case of zero population imbalance. This implies 
that we are able to find an analytical form for the pairing am¬ 
plitude at long times, although the parameters of the solution 
cannot be determined exactly from the initial conditions. 


II. MODEL 


Our starting point is the BCS Hamiltonian in the presence 
of the spin-orbit interaction in two spatial dimensions and 
Zeeman magnetic field term: 115 -■-I— I 


H = ^2 [(£k<5a/3 - h Z (J z a p) + aso(r k • 3) 

k a(3 

~ 9 E ^kt^-k.|T-k'-l- c k'T> 

kk' 


Lkcv^k^ 


( 2 . 1 ) 


where is a fermionic creation operator with momentum 
k and spin projection a, g > 0 is the pairing strength, 
r k = ( k y , —k x ), cx-so is the Rashba spin-orbit coupling con¬ 
stant, h z is a Zeeman field which determines the degree of the 
population imbalance and £k = k 2 / 2 — fi is the single particle 
energies taken relative to the chemical potential /j and we set 
the mass of the fermions to m = 1. In passing we note that 
this model, strictly speaking, is not applicable to the system 
of charged fermions since the orbital effects will dominate the 
Pauli limiting effects. 


The non-interacting part of the Hamiltonian (2.1 1 can be 
diagonalized, which yields a new spectrum 

£kA = £k — \\Jh 2 z + ( asok ) 2 , A = ±1. (2.2) 

We can now perform the unitary transformation from the orig¬ 
inal operators to new operators, which describe the fermionic 
excitations in chiral bands. The analysis of the ground state 


properties of the model (2.1 1 can be considerably simplified 
after we employ the mean-field theory approximation in the 
particle-particle channel and then make a unitary transforma¬ 
tion from the original operators c k A to a fermionic operators in 
chiral basis dkA- The resulting mean-field Hamiltonian reads: 


*=E 


kA 


ekAa^dkA - E A? 7k 0 fc«kA«- 


kA 


kA 


A 


(2.3) 


\ ^—■*, -—• 

- — 2_^ VkOka-kxa kJ + h.c. 


kA 


Here, for convenience, we introduced the following momen¬ 
tum dependent functions: ? 7 k = exp[i tan~ 1 (k y /k x )\ and 


In the next Section we introduce the model, briefly review 
its ground state properties and derive the equation of motion 
which describe the superfluid dynamics in terms of real func¬ 
tions. In Section III we analyze the possible steady states 
which appear as a result of quench in the pairing strength 
for equal atomic populations. In the first part of Section 
IV we discuss the steady state diagram for the quenches to 
the state with zero population imbalance, while in the sec¬ 
ond part present the results of the numerical simulations for 
the quenches into a state with non-zero population imbalance. 
Secton V is followed by the concluding discussion of our re¬ 
sults. Lastly in Appendix A and Appendix B we provide the 
details on the derivation of the equations of motion. 




f? k 


asok ~ _ h z 
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(2.4) 


Formally, the model (2.3 (is analogous to the model discussed 
by Sato et al. 15 ®. The crucial difference in our case, however, 
is that the pairing gap A is not proximity induced and instead 
must be determined self-consistently: 


A = ?/k A0 fc (d_kAa k A) + 0fc(a-kAa k A) 


(2.5) 


kA 
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FIG. 1: (Color online) Dependence of the pairing amplitude A, 
chemical potential g and spectral gap E gdp = _E k _ 0A __|_ (all in ar¬ 
bitrary units) as a function of the Zeeman field hz determined by 
the numerical solution of the self-consistency equations ( |2.13|2.14^ . 
Note that the superfluid becomes gapless while the order parameter 
remains finite at some critical value of the field = \/ g 2 + A 2 . 
These results correspond to the following choice of the parameters: 
n c = 0.125, ef = 2irn c = 0.785 and aso = 0.752. 


The mean-field Hamiltonian (2.3 i can be diagonalized. We 
find that the single particle spectrum consists of four bands 
u>± (k, A) = ±-EkA with the following dispersion 


Ek\ = 


& 


Rl 


A 2 — 2Af? k v/ 


02 A2 


1/2 


( 2 . 6 ) 


2(— Q k A x , — Q k A y , 0) can be interpreted as an ’’induced 
magnetization” since its xy-components vanish for hz = 0. 
Naturally, equations \2.1\ have the form of the Bloch equa¬ 
tions for the BCS superconductor when hz = 0. The first two 
components of B^ x are determined self-consistently by 


A x(t) iA y (t) — g 'y ' [0 fe S kM (t) E 

k/i 


( 2 . 8 ) 


where we have adopted the usual notation = S x ± iS y . 
Equations of motion for the components of vector T kA (£) are 

d t Ll x = - 2 e k L y kX - Q k A y (t) [S£ A + S^] 

— 20 fc A x (£)T k , 


d t Ll x =2e k L£ A + 0*A x (i) [S£ x + S^\ 

— 2Q k A y (t)T k , 

dtL^ x +2\R k T] i (t) + 0fcA x (t) 

— 0feA 2/ (f) 


(2.9) 


^ - ^A 


QX QV 

* k A - * kA 


= 0 , 


where e k = k 2 /2. Note that as it follows from these equations 
L kA = IRj and also L kA = -T‘-. Finally, the last equa¬ 
tion of motion which determines the evolution of the auxiliary 
variable T k reads: 


<9tT k + Bk X (t) ■ 4a( f) 


— 2e k L kA . 


^^7n k (f) • S kX (t) 

A 


( 2 . 10 ) 


Before we discuss the ground state properties of the model 
( |2.3| ), we first introduce the auxiliary functions which are anal¬ 
ogous to the pseudospin variables for the BCS model. 


A. Equations of motion 

In this Section we list the equations of motion (EOM) 
which will allow us to study the dynamics of the pairing am¬ 
plitude in the collisionless regime. EOM can be obtained from 
the corresponding EOM for the single particle propagators, 
which can then be cast into the form of the EOM analogous 
to the Bloch equations for the magnetic moments in exter¬ 
nal magnetic field. As a reader may have already guessed, 
there should be ten equations of motion overall: six equations 
describe the Cooper pair dynamics on each of the two chiral 
bands A = ±, while the remaining four appear as a result of 
non-zero Zeeman field. The details on the derivation of the 
equations of motion are given in the Appendix A, so here we 
provide the final results. The first six equations are compactly 
written as follows 

<9 t 5 k a = B kX (t) x S kX (t) + m k (f) x 4a(*) (2.7) 

with f? kA = 2(—0fcA x , —Q k A y , £ kA ) is an effective 

field around which S is precessing and vector m k = 


As we can immediately observe from these equations of mo¬ 
tion, in the absence of the Zeeman field the first six equa¬ 
tions decouple from the rest and become equivalent to the An¬ 
derson equations of motion for the pseudospins in the BCS 
mod cl. 1221211 Thus, based on this observation we conclude that 
the evolution of 5 k A (f) can be determined exactlyP2EH How¬ 
ever, for the general case of nonzero Zeeman field, one needs 
to resort to the numerical solution of the equations above for 
the dynamics initiated by a sudden change in the parameters 
of the model, such as pairing strength g , Zeeman field hz 
or spin-orbit coupling aso- I n what follows we specifically 
study the quenches of the coupling constant and Zeeman field. 


B. Initial conditions 


Let us write down the expressions for the auxiliary func¬ 
tions Sk X (t), L^ x {t) and T k (£) at time of a quench, t, = 0. In 
what follows we only focus on the case when the system is ini¬ 
tially in its ground state. Then, the initial momentum distribu¬ 
tion for these variables directly follows from the equations of 
motion ( 2.7|2.9|2.10| i. Without loss of generality, we assume 
that initially the superfluid order parameter is real, A x = A, 
A y = 0. Employing the relations between the single particle 
propagators, evaluated at equal times and auxilary functions 
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above, for the x components of S k x(t) and L k x(t) we find 


-Sk*(0) = 


LZ X ( 0) = 


0fcA (E kX E kX + £“ x 


A 2 ) 


2E kx E k j(Ekx + E kx ) 

0fcA[£'kA-E'i c X e kA^kA "I” A 2 ] 

2E k xE kX (E k x + -^kl) 


( 2 . 11 ) 


while iS kA (0) = L kA (0) = Tk = 0. Reader can easily check 
that in the limit hz = 0 we recover the expression for the 
Anderson pseudospin in the BCS model. Consequently, in 
the limit of aso = 0 we naturally find S^aW = 0 while 
££a(0) = <Pk&/2E k with E k = x/(e k - \i) 2 + A 2 and 
<Fk = £k[l + sign(£k - hz)]/{E k + hz + \E k — A z |). 
Similarly, for S£ A (0) and L kA (0) we 0 b ta i n: 


jz _ QfcQfcA"(f kA £kA) 

kAl J 2E k xE kx {E k x + E kx V 

^a (0) = -^^a (0) + |^a (0 ). 


( 2 . 12 ) 


One can easily check that in the limit hz = 0 we recover the 
usual expression for S kA (0) in the reduced BCS model. For 
the case of finite Zeeman field and no spin-orbit coupling L kx 
is zero, while A[S£ a ( 0)] qso=0 = -(e k - a)[££a(°)]c«so=o- 
In this case there is a similar decoupling in the equations of 
motion and we only need to solve six dynamics equations 
instead of ten. As it turns out, the pairing dynamics in this 
case can be found exactly. 47 Next, we discuss the ground state 
properties of our mean-field model. 


where we used the relation 0^ + 0 2 = 1, n c = £f/2tt is a 
particle density per spin in two dimensions and £p is the Fermi 
energy. We analyze both of these equations numerically and 
present the results of our analysis on Fig. [I] Perhaps the most 
remarkable feature of our results is the vanishing the spectral 
gap Eg a p = Ak— o.a=+ at some critical value of the Zeeman 
field hzc = + A 2 , while the pairing amplitude remains 

finite. This effect is well understood: it signals a topologi¬ 
cal phase transition at which the winding number W changes 
from W = 0 to W = 1 (for related discussion see e.g. Ref. 
Il34l and references therein). The change in the winding num¬ 
ber reflects the appearance of the Majorana gapless chiral edge 
modes in a sample with boundaries. 


III. QUENCH OF THE PAIRING STRENGTH IN THE 
MODEL WITH ZERO POPULATION IMBALANCE 

In this Section we consider the pairing dynamics following 
the sudden change of the pairing strength for equal atomic 
populations, hz = 0. In this case 0^ = 1 and 0^ = 
0. We will mainly focus of the details of the steady state 
’’phase diagram” ignoring another aspects of the problem such 
as long-time asymptote of the pairing amplitude and steady 
state quasiparticle distribution function due to the similarity 
with the corresponding problem discussed in great details by 
Yuzbashyan et al. l40l . 


A. Lax vector 


C. Ground state 

In the absence of spin-orbit coupling superconductivity be¬ 
comes energetically unfavorable when the magnitude of the 
Zeeman field is hz > a/ 2 A known as Clogston-Chandrasekar 
criterion.®® Nonzero spin-orbit coupling, however, leads to 
the mixing between singlet and triplet components in the 
anomalous Gor’kov correlation function^ and superconduc¬ 
tivity extends to much higher values of the Zeeman field. 

The value of the pairing amplitude in the ground state is 
determined from the solution of the self-consistency equation 
(|2.8|>. Taking into account equations (|2.11|i above, we find 


Here we will introduce quantities, which we will later use 
to analyze the steady state dynamics of the condensate. The 
Lax vector for our problem is defined according to: 


A(u) = Y 

kA 


S kA _ e/ 

u - £ k A 9 


(3.1) 


Equation of motion for the Lax vector follows directly from 
the equations of motion for the pseudospins S k x’- 

d t l{u) = [—2A(i) + 2 ue z ] x C(u). (3.2) 


The square of the Lax vector is conserved by the evolution 


1 = E ^ E k\ + A 2 + 0fcC kA + 0fc£kA£ k A (2 13 
9 k\ 2E kx E kJ (E k x + E kx ) 

In addition, we need to compute the value of the chemical po¬ 
tential /i in the ground state. The equation for the chemical 
potential is obtained from the standard expression for the par¬ 
ticle number in terms of the functions S kx . We find: 


2 n r = 


kA 


s k x{E k xE k - x +el^ + Qltf) 


+ 


2E k xE kx {E k x + E kx ) 


0 fcA 2 £ kA 


2E k xE kx {E k x + E kx ) 


(2.14) 


i\u) = 4 
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£ 

pA 


2H 


pA 


s: 


pA 


£pA ^pA)^ 


(3.3) 


where we have introduced 


'Hpx = 


S 


pA 


J q/^ 




- (3.4) 
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Following the arguments of Ref. go) we immediately con¬ 
clude that the dynamics governed by the mean-field Hamilto¬ 


nian (2.1 1 with hz = 0 can be determined exactly. 

Our main goal in this Section is to determine the steady 
state phase diagram, which we will plot in the plane of initial 
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and final values of the superfluid order parameters, Aoi and 
Aof, just like it has been done in earlier works. 32 34 ® 

As it has been extensively discussed in Ref. ESI, in the 
thermodynamic limit the imaginary part of the complex roots 
of the spectral polynomial determine the value of the pairing 
amplitude in a steady state. Let us compute the roots of ( |3.3| l 
for the initial configuration of the pseudospins. It follows: 


Next, we introduce an integration variable e = k 2 /2 — asok, 
so that: 


k± (e) = aso 1 ± W1 + 


2e 


l so, 


(3.13) 


q- 


£ x (u, 9i ) = V " kA = A 0i £ 0 (u), (3.5) 


For the first integral in (3.12» we need to pick k_(e) while in 
: fc+(e). It fol 

de 2 a S o 


the second integral we pick fc + (e). It follows: 
o 


E^ ek +) = 


where 


>/ 2 


2tt 


SO 


2e 


A)0) = 


(3.6) 


kA 2(M-e kA )v / (£kA -/x) 2 + A 2 

Similarly, £ y (u,gi) = 0 and 

£ z (u,gi) = -(u- g)£ 0 (u). (3.7) 

Thus, Eq. ( |3.3| ) becomes 

[(u-g) 2 + A 2 1 ]£ 2 (u) = 0. (3.8) 

Clearly, the equation ( |3.8[ > has the complex conjugates pair of 
roots: 



F(e) 


- F(e). 


(3.14) 


The contribution from the chiral band A = — is trivial and it 
yields: 

OO 


E^( £k -) = 


Thus, Eq. p.l l[i becomes 



F(e). (3.15) 


de 


u 0 ,± = M ± iA 0 


(3.9) 


and the imaginary part of iio : ± gives the value of the pairing 
amplitude We also define a spectral polynomial 

Q 2 N+ 2 (u) = g 2 n(« - £ pa) 2 • A' 2 (u), (3.10) 

pA 

where N is the total number of distinct single particle energy 
levels £pA. Since we are considering the case when the pairing 
strength changes abruptly from gi —> gf, we set g = gf in 
Eqs. > [3.1|3.10| ). 


u- g =F *A()i 
0 


2{u - e)yj(e- g) 2 + A 2 


>/ 2 


asode 

2 V a so + 2e(u - e)^/(e - g) 2 + A 2 { 


= 0. 


(3.16) 


where /3 = f3/2vp and lod is the bandwidth. Naturally, when 
oso = 0 we recover the equation for the Lax roots in the BCS 
model. Although in the subsequent analysis we can safely 
take wc -> oo, however, in numerical calculations we have to 
keep the bandwidth finite. 

We are interested in finding the values of /? for which the 


B. Roots of the spectral polynomial and steady state diagram equation (3,16) will have two pairs of complex conjugated 


For the case when the coupling is changed instantaneously, 
the complex roots of Eq. (|3.3| > or, equivalently, the roots of the 
spectral polynomial (3.10| with g = gf can be obtained from 


P 


roots. Let us introduce the following variable: 

u = g + v Aoi. (3.17) 

The imaginary roots which determine the value of the pairing 
amplitude in the steady state are determined by setting 


u — g =F *Aq 


kA 


2 (u - £kA)\/( £ kA - g) 2 + A 2 


= 0 , 


(3.11) 


where /3 = gj — g[ . To analyze Eq. (3.11 1 it is convenient 
to go from summations over momentum to the integration 
over energy by introducing the density of states v{e) = vp 
where vp = n c /eF and £p is the Fermi energy, n c is a par¬ 
ticle density per spin. We need to consider contribution from 
each chiral band separately. 

Consider A = + first with e k + = k 2 /2 — asok: 

a. oo 

, f kdk . f kdk 
F{e k+) = I —F(£ W+ ) + j F(e k+ ). (3.12) 


u —> u ± id. 

Using ( |3. 17| > we re-write (3.161 as follows: 

OO 

2/3(u±z) [' de 


(3.18) 



(3.19) 


~( a so^ 2 ^)/ 2 A 

de 


= 0 . 


+ ^^(v-eWe 2 + 1 














































Let us find the critical value of f3 when the imaginary part of 
v becomes non-zero for the first time. We have 
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2/3 7rt?(uAoi + n) 

V 2 + 1 ^ y/v 2 + 1 

a 2 so uA 0 i — fi)'d(2vA, +2/J + a 2 so ) 


2Aoi 


+/ 




= 0, 


de 


+ 1 J-fi/ A 0i (u - e)v/e 2 + 1 


/T2 I-- 11/A oi 


de 


^(e) 


= 0. 


(3.20) 



where we introduced for brevity function 


FIG. 2: (Color online) Steady state diagram for the case hz = 0. 
In the region I A(t — » 00 ) = 0. In the Region II A(t — » 00 ) = 
Aoq. Lastly, in the region III A(t) varies periodically with time. 
Parameters: n c = 0.125, aso = 0.75. In the limit of zero spin- 
orbit coupling the read line inside the Region III is absent. 


Let us analyze the first equation in ( 3.20| i. Depending on the 
value of v, there are two possible solutions. First solution 
corresponding to the usual BCS case: 


m = ^Vv 2 + 1, v>~n/ A,)i, (3.21) 


while v is found by solving 

7rsign(/3 e ) f°° de 

Vv 2 + 1 J-fj ,/Am {v ~ e 2 + 1 


+ 


1 a 


so 


2 Ao 


— ^/ ^Oi 

/ 


de 

W) 


= 0 


(3.22) 


still for v > —fi/ Aoi- There is, however, another solution for 
Pc given by 


\Pc\ 


iraso Vv 2 + 1 

2 y / 2Aoit7q^r| 0 ^F2/r 


-M/A oi- 


(3.23) 


IV. QUENCH OF THE POPULATION IMBALANCE 


As we have seen already, non-zero Zeeman field breaks 
integrability. Thus, for the quenches of the Zeeman field 
hz, —> hzf one needs to resort to the numerical analysis of 
the equations of motion. The main interest in studying this 
particular type of quench is mainly motivated by the existence 
of the topological transition. The task of analyzing steady 
state diagram for an arbitrary values of hz / has been recently 
accomplished by Dong et ah® However, as it became clear 
from our discussion above, for the spacial quenches such that 
hzf = 0 the problem can be analyzed analytically using the 
same method of Lax vector construction. The only difference 
with the previous analysis is that an initial pseudospin distri¬ 
bution explicitly depends on hzi- 


The value of v in this case will be given by 


nasosigniPc) 


de 


Vv 2 + ly/2AoiU + a 2 so -F 2/r 


-/V A oi 


(v - e)Ve 2 + 1 


+ • 


x SO 


2A 0i j~_ 


— li/Aa 




de 


=0 


(3.24) 


for v < —/.//A 0 ;. The results of our analysis of equations 
for the critical /3 above are shown on Fig. [2] The presence 
of the spin-orbit coupling leads to the appearance of the re¬ 
gion where pairing amplitude goes to a constant (Region II) is 
realized inside the region where the pairing amplitude period¬ 
ically varies with time (Region III). 


A. Integrable dynamics: h Z f = 0 


We start with the analysis of the expression for the Lax 
vector (3.1 1 . The expression for C z {u) can be considerably 
simplified if we take into account the self-consistency equa¬ 
tion ( 2.1 3| >. However, one needs to be careful, since at large 
fields self-consistency equation does not have a solution and 
we have to set A = Aoi = 0 in & Therefore, we have to 
consider two cases: in the first case Ao, in the initial state is 
nonzero, while in the second one hz is large enough so that 
A 0 i = 0. 


We first analyze the roots for the case of finite A U; . The 
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FIG. 3: (Color online) Imaginary parts of the roots of C 2 {u ) = 
0 and superfluid order parameter Ao; in the initial state plotted as 
a function of initial value of the Zeeman field hzi = hz for the 
quenches hzi —> hzf = 0. Note that the imaginary parts of both 
roots essentially coincide with each other for the initial conditions 
with Aoi —► 0. At hz = h c2 the second complex root appears. 
Thus, h c 2 separates the steady states with constant and periodically 
oscillating superfluid order parameter. These results correspond to 
the following choice of the parameters: n c = 0.125, Ef = 0.785 
and aso = 0.752. 



FIG. 4: (Color online) The dependence of the critical Zeeman 
field h c 2 , which separates steady state with constant and period¬ 
ically oscillating pairing amplitude, on the strength of the spin- 
orbit coupling. Note that for the case of weak spin-orbit coupling, 
the pairing amplitude will always go to a constant at long times, 
A(f —¥ oo) = Aqo. The inset shows the dependence of the real 
part of the second root on aso- These results correspond to the fol¬ 
lowing choice of the parameters: n c = 0.125, ef = 0.785. 


roots are the computed numerically from 


E 

kA 


(u — /r + *0fc A)(Ek\E kX + A" + e^) 


2 (u — /J, — £k\)EkxE k j(E\c\ + E k j) 


+E 

kA 

= E 


0 I A2 ( £ kA - £ kA) 


2 (u — H — £kA)^kA^ k y(£ ; kA + E k j) 


(4.1) 


h 2 z 


kA 


Ek\E kx (Ek\ + E kx ) 


which follows from ( |3. 1| > and the self-consistency equation 
( 2.1 3| ). We analyze this equation numerically and plot out re¬ 
sults on Fig. [3] As expected, for relatively small values of 
hz-, = hz there is only one complex root, which means that 
the steady state order parameter asymptotes to a constant. As 
the value of the field is increased further, it reaches h c2 where 
the second complex conjugated root appears. For quenches 
of Zeeman field with hz > h c 2 pairing amplitude periodi¬ 
cally oscillates in time. Our results confirm those found from 
the numerical simulations.® Indeed, on Fig. [5]we show A(t) 
found by numerically solving the equations of motions for 
various values of hz and it is clearly in agreement with our 
analysis of the Lax roots. Fig. [3] 

In Fig. [4]we also plot the dependence of h c2 on aso 1 which 
we determine by setting u = uq + iS in ( |4. 1[ ) and solving them 
together with Eqs. ( |2. 1 3|2. 14| >. As one may have expected, 
h c 2 ex asoPF- Furthermore, the fact that we do not find a so¬ 
lution for small aso is i n qualitative agreement with an obser¬ 
vation that the steady state with oscillating pairing amplitude 
generally appears for moderate to strong quenches. 

Next, we would like to show that no more complex roots ap¬ 
pear at large fields when Ao, is infinitesimally small. First, let 
consider the case when the self-consistency equation (|2.13 1 


does not have a solution and, as before, we set u = uq + i5. 
Then, in the equation for the Lax roots C z (u) = ±iC x (u) we 
can consider the real and imaginary parts separately. Equation 
for the imaginary part is satisfied only if uq < — hz, while the 
equation of the real part reads 


- + yp( sign(a ‘ A) 1 = 0 , 

9 “ \Uq — p — £kA, 


(4.2) 


where P stands for the principal value. We have analyzed 
this equation and did not find a value of coupling g consis¬ 
tent with the zero value of the superfluid gap. We reach the 
same conclusion from the analysis of Eq. ( |4. 1 [ > for the case 
when An, is small enough so it can be neglected. To summa¬ 
rize, we find that for the quenches of the Zeeman field from 
some finite value hzi to zero, there are only two steady states 
possible at long times: in the first one A(t) asymptotes to a 
constant, while in the second one A(f) continues to oscillate 
periodically. 


B. Analytical solution for the pairing amplitude 

In this subsection we derive the analytic expressions for the 
pairing amplitude A(t) in a steady state. Our discussion here 
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follows closely the related discussion in Refs. 021401 . 


The Lax vector describing the reduced solution is 



FIG. 5: (Color online) Results of the numerical solution of the equa¬ 
tions of motion for the A (t) in the exactly integrable case when 
hzf = 0. These results correspond to the following choice of the 
parameters: n c = 0.125, ef = 0.785 and h c 2 = 1.02£f- 


^red(^) — 

j 

1 + E 
\ P A 

f m _> 

£m{u) = 

a i 

^ u — 

U= 1 




(4.3) 


Here C m {u) is a Lax vector for a reduced system, time- 
dependent vectors oj and parameters d p x,£j need to be de¬ 
termined. As it can be easily seen, vectors 5j satisfy the same 
equations of motion as original pseudospins S p \. The param¬ 
eters of the reduced Lax vector are chosen such that 


£(u) = C reA (u). (4.4) 


Therefore, the equation of motion for vector £ le d(u) is the 
same as the one for C(u): 


d t C red (u) = [—2A(f) + 2 ue z \ x £ red (u), (4.5) 

where we use A = (A x ,A y ) for brevity. By matching the 
residues at u = e p \ and at u = ej we find the following set of 
relations: 


E 



= -1, 


d p X Am (^pA) — ^pA- 


j = 1, ..., TO, 


(4.6) 


In the thermodynamic limit it is possible to find the reduced 
solutions that have the same integrals of motion as the so¬ 
lutions for the quench ed dy namics, i.e. they have the same 
C 2 (u). Thus, equation (4.3 1 becomes 


i+Y, ( " x{ex)dx{ex)dex = 


u-e \ 


I C 2 (u) 


(4.7) 


where £(u) = ±1, vx{t) is the density of states for the chiral 
band A and C 2 [u) is determined by the initial conditions. By 
setting u = e ± iS we can immediately determine dy(e) 


<(e) ( V£ 2 (e~) _ V^ 2 (£+) \ 

2nv x {e) {^J^j J 


(4.8) 


with e± = £ ± id. In what follows, we will derive the explicit 
expressions to determine parameters dj and ij :l (j > 1), which 
define C m (u), in terms of the complex roots of C 2 {u). 

a. m = 1 solution. This is the case of the one-spin so¬ 
lution. The expression for L m -\ reads: 


Steady states with constant and periodically oscillating 
pairing amplitude can be described analytically by construct¬ 
ing the Lax vector for an effective m-pseudospin system. The 
Lax reduction procedure states that at long times the dynamics 
of a superfluid is governed by a dynamics of only few gener¬ 
alized pseudospin variables, which we will denote by dj. E 2 ED 


C m =i(u) = —-— • (4.9) 

u-e 1 g 

The relation between <j\ and A follows directly from the self- 
consistency equation ( |2.8| > and Eq. ( |4.6[ > : 

Alt) = gdxlt), 


(4.10) 
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which also implies that of remains constant. Using the equa¬ 
tion of motion for the Lax vector ( |4.5[ > together with ( |4. 1 ()[ > we 
can now solve for A(f): 

A(f) = A oc e 2i,J ” x,t ~ ilpo , (4.11) 

with Hoo = ei + gof and ipo is an integration constant. Pa¬ 
rameters {Aoqj/lIoo} can be expressed in terms of the roots 
for the square of the reduced Lax vector. Recall, that in the 

thermodynamic limit these roots are the same as the roots of 

-.2 

£ (it) = 0 by construction. As we have seen in the previous 
section, when hzi < h c 2 there is only one pair of complex 
conjugated roots, which we denote u± = uu ± iuu. Taking 
the square of the both parts in Eq. ( |4.9[ > and regrouping the 
terms in the right-hand-side yields: 

tilt — A^oo 1 — Aqq. (4.12) 


follows that £ is conserved provided the terms containing 
A(i) drop out from (4.15 1 . In turn, this is only possible for 
erf 2 cx a 2 . Therefore, we write 122LO] 


cr~ = + 61 , erf = — n 2 + b 2 , 

9 “ 9 

where coefficients 01,2 and 612 satisfy 

eiai + €20.2 = 2(ei6i + €262) = £, 
Ol = —02, bi + b 2 = cr. 


(4.16) 


(4.17) 


Importantly, by a virtue of the second equation (4.6 1 we obtain 
the following ansatz for the original variables: 


Sp\ — a pA ^ 2 + b p \. 


(4.18) 


Thus, in agreement with earlier results we find, that the imag¬ 
inary root of C 2 (u) = 0 determines the value of the pairing 
amplitude at long times. 



FIG. 6 : (Color online) Dependence of the roots ei, e 2 and e 3 of the 
qubic polynomial Ps(w), Eq. 1 4.25 1 on the value of the imbalance 
h z . n c = 0.125, e F = 0.785. 


Furthermore, equatio ns of motion for the two remaining com¬ 
ponents of S p \ - Eq. (2.7 1 with m = 0 and 0^ = 1 - yields: 


S~^-S+, e"** = 2 *o pA fi ) 


pA 




(4.19) 


where we use the notation Sp X = S* x ±iSp X . After a series of 
algebraic manipulations identical to the ones in Refs. 1321401 . 
we find the following equation for fl: 


n 2 + rr + 


26pA 

^pA 


+ 4£pv ) ^7“ — 44ep>,f2 


A 2 


blx-S 2 pX 


(4.20) 


= 0, 


where A is a function of (1 given by 

A = 2 ha£1 + 


(4.21) 


and ha, ka are arbitrary real constants. Since the same equa¬ 
tion for fl is found by considering the equations of motion 
for the variables < 71,2 we conclude that the coefficients in Eq. 
(4.20 1 must be independent of p and A: 


b. m = 2 solution. This is the case of the two-spin so¬ 
lution with the reduced Lax vector of the form 


£'m= 2 (^) 

| (J 2 e 2 

(4.13) 

u-e 1 u e 2 g 


and 



A = g ■ (<xi + o 2 ), 

A x (t) - iA y (t) = 

(4.14) 


where in the second expression ( .l = A and $ is the phase 
of the pairing amplitude. The dynamics of the variables t?i ,2 
is governed by the following two-spin Hamiltonian: 


~~—f 2(e p x — ha) 2 — 2 p, 

tlpA 

h 2 - Q2 

°pA ^pA 


X pA 


(e p a - ha) = 4\. 


Thus the differential equation for f l(t) becomes 

Cl 2 + n 4 + 4 pU 2 + § + 4 * = 0. 

Solution of this equation is® 


(4.22) 


(4.23) 


Hm =2 = 2{e 1 o z 1 +e 2 o^)-2A-o. (4.15) fi = s/A 2 +e u A = A + dn[A + (f - t 0 ),k'}, (4.24) 


The ^-component of S = <T \ + a 2 is conserved by evolu¬ 
tion governed by the reduced Hamiltonian (4.151, which re¬ 
flects the total particle conservation. In addition, the total en¬ 
ergy £ must be conserved by the evolution. Given the self- 
consistency condition (4. 14[> for the reduced Hamiltonian, it 


where dn is the Jacobi elliptic function, k' = A_/A + , A 2 _ = 
e 2 — ei, A+ = e 3 — ei and the parameters ei 2,3 are the real 
roots of the qubic polynomial 

Pa(w) = w 3 + Apw 2 + A\w + k\. (4.25) 
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The last step is to match the coefficients in the polynomial 
( |4.25| ) with the values of the complex conjugated roots ap¬ 
pearing for hzi > h c 2 , Fig. |3| To do that, we will employ 
the relation (4.6 1 . First we solve Eqs. (4.22 1 for a p \, b p \. We 
find 


5, 


^pA 


b p x — 


pA 


2 \/[( e pA - Ha) 2 - P ] 2 - K A (e pA - ha)-X 

[( e pA — Ha) — Pl^pA 


\/[( e pA - Ha) 2 - P ] 2 - ka(£ p a - Ha) - X 


(4.26) 


Similarly, the coefficients 01 , 2 , 61,2 of the reduced solution 
(4.16) are found using the conservation laws (4.17 1 : 


ai,2 = ± 


2(ei — £ 2 )’ 


& 1,2 = ± 


£ - e 2 ,icr 
£i - £2 


(4.27) 


Using these expressions, let us mat ch the pre-fact ors in front 
of f l 2 after we use Eqs. (4.16 4.18 1 together with ( 4.26|4.27 1 
in second equation in (14.6 1 for the z-components of C m and 
S p \. We find: 


^pA 


(^pA ^l)(^-pA ^2)^pA 


9 \Z[(s P A - Ha) 2 - P ] 2 - ka(£ p a ~Ha)~X 

(4.28) 


On the other hand 


d p \ — 


5, 


pA 


\/£'rn=2 ( e pA ) 


(4.29) 


Introducing the spectral polynomial Q<i{u) similar to (3.10»: 

^2 


Qa{u) = g-(u - £i) 2 (u - e 2 ) 2 • C m=2 {u). 


(4.30) 


If we now compare ( |4.30[ ) with (4.28 1 we immediately identify 
Qa{u) with 


Qi{u) = [(it - ha) 2 - p ] 2 - ka{u - ha) - X- (4-31) 


Furthermore, since in the thermodynamic limit the complex 
roots of Qi{u) must match the complex roots of C 2 {u), we 
can express all the parameters ( |4.31[ ) in terms of two pairs of 
complex conjugated roots «i 2 = iti, 2 t + *iti, 2 > : 


HA = 


itit 


It 2t 


p = 3 h 2 a ~ 2u lt u 2x - 


U lt + U ii + U 2t + »2i 


KA — 2lti r (lt 2 r + U 2 i) 4” 2lt 2r (lti r + Itji) 

+ 4 ha(p-Ha), 

X = KAHA + (ftA - P) 2 
- (u 2 r + u\ { )[u\ x + u| t ). 


(4.32) 


We plot the dependence of the roots of Ps{w) (4.25 1 on Fig. 
[ 6 ] Note that e\, e 2 and e 3 are small for hz ~ h c2 . It was noted 
in Ref. [40]for the quenches of the detuning frequency across 



FIG. 7: (Color online) Results of the numerical solution of the equa¬ 
tions of motion for the A(t) in the general, i.e. non-integrable, case 
hzf 7 ^ 0: (a) hzj = 0.9ef; (b) hzj = 0.5ef; (c) hzj — 0.25sf 
and (d) hzf = O.Ief- The values of the remaining parameters are: 

hzi = 1.85e F , n c = 0.125, e F = 0.785. 



FIG. 8: (Color online) Same as Fig. [7] with (a) hzf = 1.25 e F \ (b) 
hzf = 1.15ef; (c) hzf = 1.1 ef and (d) hzf = 0.95ef. 


the Feshbach resonance, the value of ei serves as a measure 
of the deviation from the weak coupling limit when |ei| <C 1. 


To summarize, the equations ( 4.32) together with ( |4.23 1 
provide exact description of the order parameter dynamics in a 
steady state determined by the two pairs of the complex con¬ 
jugated roots of the spectral polynomial. In particular, the 
pairing amplitude is given by 


\A(t)\ = \Je\ + A“_dn 2 [A + (f - f 0 ), k'} , (4.33) 

where the parameters entering into this expression are given 
above, ( |4.24| >. Note that parameter e\ is close to zero only 
when hz ~ h c2 . It is somewhat surprising to find that |A(t)| 
is described by the weak-coupling solution—! 

|A(f)| oc dn[A+(f — to), k'] 


(4.34) 
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only at lower fields. 


C. Pairing amplitude dynamics with finite population 
imbalance 


Here we will discuss the dynamics initiated by the quenches 
of the Zeeman field, so that hzf 7 ^ 0. Since the dynamics 
governed by the Hamiltonian is non-integrable, we have 
to resort to the numerical analysis of the equations of motion 
( 2.7|2.9|2.10| i. Our main motivation for this part was to check 
whether the steady state with the periodically oscillating pair¬ 
ing amplitude also extends into a non-integrable region of the 
parameter space. 

The time evolution of the pairing amplitude following the 
quench is shown in Figs. 0 and [ 8 ] We see that for certain 
values of Kz/ef the order parameter magnitude A(t) shows 
oscillations with several frequencies and its amplitude is not 
constant at long times (at least up to the longest time scales we 
were able to achieve with our numerics). However, note the 
striking difference between the dynamics in Fig.[7]and Fig. [ 8 ] 
when hzf exceeds the value of /i c3 « 1 . 02 £p provided hzi = 
1.85ef. the pairing amplitude shows regular oscillations with 
constant amplitude. This behavior is characteristic of A(f) 
which is found in exactly solvable limit. 

To get further insight into the origin of this behavior, on 
Figs. [9] and 10 we plot the single particle energy depen¬ 
dence of the auxiliary functions L(e,i) and T(e. t) at long 
times when hzi < h c 3 and hzi > h c 3. For these plots the 
regular oscillatory behavior of A(f) becomes clear since for 


hzi > h c 3 equations of motion for the functions Sk\(t) de¬ 
couple from the remaining four equations of motion (|2.9|> and 
CT . Lastly we make one more observation: this dynamical 
decoupling happens exactly when the system goes through the 
Floquet topological transition^lSI corresponding to the transi¬ 
tion from topologically trivial Floquet spectrum to a steady 
state with topologically non-trivial Floquet spectrum. How¬ 
ever, the detailed analysis of this transition goes beyond the 
scope of this paper and we leave it for the future publication. 


V. CONCLUSIONS 

In this paper we have analyzed the far-from-equilibrium 
pairing dynamics of the spin-orbit coupled fermions in 2 d with 
population imbalance. Specifically, we have considered two 
separate cases. In the first case the dynamics is initiated by 
the a sudden change of the pairing strength. We found that 
the steady state with periodically varying pairing amplitude is 
realized in much narrow regions of the steady-state phase di¬ 
agram compared to what happens when spin-orbit coupling is 
zero. 

Exact integrability of the problem with zero imbalance im¬ 
plies that we can also provide analytical description for the 
dynamics initiated by a sudden change of imbalance to zero. 
We find that when the initial value of the imbalance field hz 
exceeds some critical value /i c2 , the steady state with period¬ 
ically oscillating pairing amplitude is realized and determine 



FIG. 9: (Color online) Energy dependence of L(e) and T(e) at 
time t5 = 2.05 ((5 is a level spacing). The values of the remaining 
parameters are: hzi = 1.85ef, hzf = 0.5ef, n c = 0.125, A = 
10ef and ef = 0.785. 



FIG. 10: (Color online) Same as Fig. [9] with hzf = 1.25 ef. In 
contrast with the Fig. Mwe see that all components of L(e, t 00 ) 
as well as T(e,t —> 00 ) are vanishingly small for all single particle 
energies. 


an analytical expression for A(t). 

Perhaps our most interesting result is our finding of the dy¬ 
namical decoupling for the quenches to finite values of the 
population imbalance. Specifically, our numerical analysis 
of the equations of motion showed that when final value of 
the population imbalance exceeds some value h t , the pairing 
amplitude is determined by a reduced number of the ’’pseu¬ 
dospin” variables. Interestingly, the value of the h t is a crit¬ 
ical value separating the regions of topologically trivial and 
topologically non-trivial Floquet spectrum . 34 Whether topol¬ 
ogy plays a defining role in the above mentioned reduction or 
it is just a mere coincidence is an exciting issue, which we 
leave for the future studies. 
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Appendix A: Equations of motion for the single particle 
correlators. 


In this section we will analyze the ground state properties 
of the Hamiltonian ( |2.3[ ) using the equations of motion of the 
single particle correlators. The main idea is to derive the set of 
the self-consistent equations describing the collisionless evo¬ 
lution of the pairing amplitude. 

Consider the equations of motion for the fermionic opera¬ 
tors. We have 


= EkAOkA - 




-pAj ’ 


/ r^^kA^) — ^P^^kA F ? /kA AOfcU—kA T" ©A;U_ kA 


(Al) 


Next we introduce the following correlation functions which 
are diagonal in new basis: 

G P x(ti,t 2 ) = -i (f (a P A0i)Op A {>2)) ) , 

F p \(t ll t 2 ) = -i\v P (t {a-px{t\)a-px{t 2 ))) , 

' (A2) 

Gpxih, t 2 ) = -i (f (aLp A (/i)a_pA(/ 2 )) ) , 

F p x{ti,t 2 ) = —iXr/p (f (aLpA^OapA^))) 


Similarly, we introduce the ”off-diagonal” correlators which 
account for the scattering of fermions between the two chiral 
bands: 


r pA (h,t 2 ) = -iX (f (a pX (ii)a^(i 2 )) ^ , 
® p x(G,t 2 ) = —irjp (f (a P A(fi)a_ pA (f 2 ))) , 

f P A(fi, t 2 ) = ~iX (f ^d^ p -(fi)a_ P A(f 2 )) ^ , 
$ p a (h,t 2 ) = -ii 7 * (f (a^ pX (ii)^ A (f 2 )) ^ . 


As a next step one can derive the equations of motion for these 
correlation functions using (All. 


c. Equations of motion for the diagonal chiral correlators 
For the diagonal in A correlation functions above we have 




Gkx{ti,t 2 ) 


+ A 


Bfc-^kA^l) f 2 ) + Bfc'&kA^l, f 2 ) 


<5(fl — t 2 ), 



r k A(G,/ 2 ) 


-A 


Bfc$kA(fl,f 2 ) — Bfc-FkA^li i 2 ) 


= 0 , 


(A4) 


Similarly, the equations of motion for the anomalous correla¬ 
tion functions (|A3|> are: 


-f EkA^ F kx {ti,t 2 ) 

QkGux(ti,t 2 ) + B fc r k A(/i,/ 2 ) 


+ A 


l ff ti + e kA ) $ kA(G,/ 2 ) 


- A 


BfcFkA(ii,i 2 ) — BfcGkA(/i,/ 2 ) 


* Qt 1 ~*~ £ kA ) r kA(fl,/ 2 ) 


— A 0fc4>kA(/i, f 2 ) — QkFkx(ti, t 2 ) 


= 0 , 


= 0 , 


= 0 . 


(A5) 


In equilibrium, all these correlation functions depend on [t\ — 
t 2 ) only, so we can perform the Fourier transform and com¬ 
pute them explicitly. It follows: 


TkA(w) — 


(tU + £ kA )A@fcf ? kA(w) — A 2 0fc@A-GkA(u-’) 


w2 - e kA “ A 2 ©! 


^kA(w) = — - 


(w — £ kA )A0 fc G k A(cu) + A 2 0 fe 0 fc F k A(w) 


w2 ^ £ kA^ A20 fc 


(A 6 ) 


where we assume that A = A, i.e. in the ground state the 
pairing amplitude is real. We have 


G kxM 


-PkA(w) 

r k A(w) 

^kA(w) 

fkA(w) 


( W + £kA)[0^A 2 +£ 2 x -u; 2 ] 
El x )(^ - El-) 
0 2 A 2 (w + £ kX ) 
(^ 2 -^a)(^-^a)’ 

©feA^+^-cu 2 ] 

(“ 2 -EL)y--Ei x y 

2XQ k Q k A 2 R k 
(w 2 — f? kA )(w 2 — Ely) ’ 
0fcA[A 2 + (£ kA — ut)(u> + EkA)] 

2xe k e k A 2 R k 
~ (w 2 — El x )(u} 2 — G 2 -) 


(A7) 
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The last expression follows from the symmetry properties of 
the corresponding equations of motion. Note that from these 
expressions it follows 


= $ kX (-w), r kA (w) = -r kI (w), 

f kA (w) = r kX (w). 


(A8) 


To compute the averages which enter into the self- 
consistency equation which determines A, we employthe 
Matsubara frequency representation w — t iu> n . Then per¬ 
forming the summations over the Matsubara frequencies and 
take the limit T —>■ 0. The resulting functions of momen¬ 
tum are listed in the main text, Eqs. (|2.1 1|2.12). Note that 
L 




oc l a kA a kx) is generated already within the mean-field 
theory despite the fact that the terms proportional to aL-ak A 
do not enter into the Hamiltonian. In what follows, we will 
also consider function 


?k — [r kA (* w ii) + rkA(*CU ra )] 


(A9) 


which is zero in the ground state, however it is generated dur¬ 
ing the evolution. 

Our goal now is to derive the equations of motion for all the 
correlation functions above as a function of 


t = 


t\ + f 2 


(A10) 


Let us introduce the following functions 


SZx{t) = g [G kA (f)-G kA (f) 

Su X (t) = S£ x (t) - iSl x (t) = -iF kX (t), 


(A14) 


and 


SkA(*)=SkA (t) + iS y kx (t) = -iF wx (t) 


— 2 r kA (f) - r kA (f) 

K a(*) = ^kA(t) ~ iLt x (i) = -^ kA (f), 

Li x {t) = Ll x {i)+iLl x {t) = -i^ x {t) 


(A15) 


In terms of these new functions, self-consistency equation for 
the pairing amplitude reads 


A 0) = 9 [ 0 fc' S 'k M ( < ) + Q kL^(t) 

k n 


(A16) 


d. Equations of motion for the off-diagonal chiral correla¬ 
tors The remaining equations of motion for the components 
of L can be derived in the same way. Let us obtain the equa¬ 
tions of motion for T kA (f). In what follows the only relation I 
will use is 


Since both normal and anomalous correlators ( |A2|A3| > depend 
on t = t, \ t‘ 2 ~ but the order parameter A(t) is a function of 
total time t only. Thus, in what follows we consider r = 0. 


Lrom equations of motion for the fermionic operators (All 
and ( |A2|A3| > it follows 


i—r.G\t X {t) + A(f) QkFk X (t) + Qk®k\(t) 

at L 


- A(f) 0fci 7 k A(f) + 

d 


i-jj. — 2e kA J F kA (f) + A (t) |©fc 


= 0, 

G kA (f)-G kA (f) 


+©fc 


r kA (f) — r kA (*)]} = o. 


(All) 


Similarly for the remaining two correlation functions which 
are diagonal in new basis 1 find 


F kA (i) — FkAtf)) — $ kA (i), 

G kA = -G, 


(A17) 


kA- 


The validity of these relations will be proven when we analyze 
equilibrium. We need to keep in mind, however, that given the 
second relation we expect that equations of motion for L kA (t) 
should not depend on chiral band index A. The equations of 
motion for the correlator < E > kA (ii, is) are 


^+ e kx) $kA(tl,t 2 ) 


- A 


©fcr k A(*i,*2) — e fc G kA (fi,f 2 ) 


+ £kA^ d*kA (tl 5 1-2 ) 


— A 


©ftTkA^l; tf) + 0fcG kA (fl, f 2 ) 


= 0 , 


= 0. 


(A18) 


i^G kA (i) + A(f) OkFk\{t) + ©fc$ kA (f) 


- A(f) 0fcF kA (f) 

d 


Qk&k\{t) 


= 0, 


t ^ + 2e kA j Fk X (t) + A (t) |0 fc G kA (f) — G kA (f) 


+0fc 


r k A(f)-r kA (f)]} = o. 


(A12) 


Lrom the equations for the normal propagators it follows 

G kA (f) = —G kA (t). (A13) 


where we used rj- p = - f/ p . Adding these two equations 
yields 


d 




+ 2e k J 4> kA (f) — 0feA(f) G kA (f) - G kA (f) 


= 0 
(A19) 

Lrom this equation we can immediately obtain the equations 


of motion for L kA using Eqs. (|A14 A15 1 . 


Lastly, we derive the equation of motion for 


Tic (t) = 


r k A (t) + r kA (t) 


2 


(A20) 
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Before I write down this equation, let me first obtain the equa¬ 
tions of motion for 1 kA and I\a- We have: 


' l TTT - e kA ) r kA(ti,i2) 


dt\ 


-A 


©fc^kA^li^) — ®kFkxitlih) 


-F £kA ) TkA(il, ^ 2 ) 


-A 
c 

dti 


Qk^kxit-i, t.2) — QfcF k ^(ii, i 2 ) 


+ £ kA ) r kA(il,i 2 ) 


- A 


©fc^kA^l, * 2 ) — (~h-F-k\(t-i , t‘ 2 ) 


- e kA ) rkA(fl, h) 


- A 


0fc < f ) kA(fl, ^ 2 ) — ®kF k j(t ll t2) 


= 0, 


= 0 , 


= 0 , 


= 0. 


(A21) 


where we have employed ( A17j», ry_ p = —r/ p and AA = — 1. 
Adding the first and the second equations and then the third 
and the fourth one and setting r = 1 1 — f 2 = 0 yields 


— 2 A Rk^j r kA (f) — Bfc [A4> k A(f) + A$kA(f)] 

+ 0fc [AF k A(*)+AF kX (f)] =0, 

+ 2XR k ^j f k A(f) — 0fc [A4> k A(f) + A$ k A(f)] 

+ 0fc [AF kA (i)+AF kX (t)] =0. 

(A22) 


when either aso = 0 or hz = O.Note that both L z and 7’k do 
not enter into the Hamiltonian and are generated in the course 
of dynamics. 


Appendix B: general relations between the components 
auxiliary functions in equilibrium 


We assume that in equilibrium A x = A and A y = 0. This 
implies that both at t = 0 both ,S kA = 0 and L kx = 0 in 
accordance with the self-consistency conditions. This guar¬ 
antees that seven out of ten equations ( |2.7|2.9|2.10| > for the 
components of vectors S, L and Tk are identically zero. The 
remaining three equations are 

EkASL + A • (0 fc S£A + 0 -^a) = 0, 

2e k Ll x + Q k A[SZ x + S^] =0, (Bl) 

2A R k Ll x + A {20 fc L^A - ©fc [SI a + S^] } = 0. 


Lastly, let us verify if expressions for the spin components 
satisfy (Bl 1 . For the first two equations we find: 


A 0fcSkA + ©fcTkA — ^SkA'S'] 
Q k A [S£ A + S kX ] = —2e k LkA- 
Lastly, let us check the third equation <E} : 

2Q k AL kX - O k A [S£ A + ,S- A 
40 fc 0 fc A 2 3 4 5 i7 2 

2^kA^kA( £ 'kA + E kx ) 


kA> 


(B2) 


(B3) 


where R k = sfV} 
that given the property (|A17[) we have 


al^k 2 . 


From these equations we see On the other hand 


kA 




iL 


kA- 


(A23) 


It is now straightforward to verify that the equations of mo¬ 
tion for these objects are the same as the ones listed in the 
main text, Eqs. ( 2.7||2.9 2. 10| . Thus we have ten equations 
of motion. These equations are decoupled into six plus four 


2XR k L kx — 2A R k 


0fc0fcA 2 (e kA — £kx) 
2£ ; kA£ ; kA(- E kA + E^j) 

40 fc 0 fc A 2 7? 2 


2£kA-E kA (i7kA + E kx ) 

Thus the third equation in (|B1|> holds. 


(B4) 
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